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Successive pairs of pseudo-random numbers generated by standard linear congruential 
transformations display ordered patterns of parallel lines. We study the “ordered” and 
“chaotic” distribution of such pairs by solving the eigenvalue problem for two-dimensional 
modular transformations over integers. We conjecture that the optimal uniformity for 
pair distribution is obtained when the slope of linear modular eigenspaces takes the 
value n op t = inaxi til (pf \/p — 1 ), where p is a prime number. We then propose a new 
generator of pairs of independent pseudo-random numbers, which realizes an optimal 
uniform distribution (in the “statistical” sense) of points on the unit square (0,1] X (0,1]. 
The method can be easily generalized to the generation of uples of random numbers 
(with k > 2). 

Keywords: Modular transformations. Order and chaos. Pseudo-random number gener¬ 
ators. 

1. Introduction 

Consider the standard linear congruential pseudo-random number generator 
(see ReflU for a review on pseudo-random numbers generation) 


x t +i = gx t mod p, (1) 

with p= 2 31 — 1 and g = 7 5 = 16807 (also called in this case “minimal stan¬ 
dard” □). It is well known that successive pairs (xt,xt+ i), normalized in the 
interval (0,1], show a pattern of parallel lines on the square S = (0,1] x (0,1] 13 . 
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This observation is even more general, since this pattern is observed for all 
pairs (x t ,xt+i), with l = lH However, the number of parallel lines 

over which the pairs are distributed strongly depends on the value of l. For 
instance, for l = (p — l)/2, the pairs are distributed on a single line, while for 
l = (j>— l)/2 ± 1 the pairs are distributed over so many parallel lines that the 
pairs could be considered as randomly filling the square S. This curious behav¬ 
ior has been called “order” (pairs distributed over a few lines) to “chaos” (pairs 
distributed over several lines) transition in Ref. □. This remark is of course very 
important to establish the properties of those correlations in pseudo-random 
number generation which are due to the “lattice” structure created by the prop¬ 
erties of the integer field Z. p = {0,... ,p — 1} over which map (jl|) is defined. 

In this paper we analyse the pattern of parallel lines using the theory of 
continuous bijective transformations defined by unimodular 2x2 matrices i (also 
called in mathematical papers Continuous Automorphisms of the Torus, CAT), 
which we recall in Section 2. We give a simple explanation of why the pairs are 
distributed either “orderly” or “chaotically”, by studying the dependence on l of 
the slope of Linear Modular Eigenspaces (LME), after solving the eigenvalue 
problem for CATs (see Sections 2 and 3). 

By introducing appropriate probes of the uniformity of pair distributions 
over the square S (Section 4), we are able to guess the value of the slope of the 
LME (and correspondingly the value of l) which gives the most uniform pair 
distribution. 

This study suggests the use of linear congruential pseudo-random number 
generators of the kind (jj|) to generate pairs of independent pseudo-random num¬ 
bers, by appropriately choosing the initial pair over the most uniform LME. We 
exploit this idea in Section 5, but we discover that the requirement of optimal 
uniformity in the “geometrical” sense does not produce the best pair generator 
in the “statistical” sense. Hence, we modify the pair generator by introduc¬ 
ing a variation of the LME slope, which produces a better randomization of 
the successive pairs. We compare the uniformity of the generated pairs with 
those obtained with a commonly used pseudo-random number generator ( ranO 
in Ref. i). 

In Section 6 we draw some conclusions. 


2. Modular transformations 

Continuous bijective transformations (isomorphisms) of the torus (the square 
S = (0,1] x (0,1] with the opposite sides identified) are maps defined by uni¬ 
modular (unitary determinant) 2x2 matrices A of integers B, 


(x(t + 1)\ 
\y(t + i)J 




mod 1 , 


( 2 ) 


where x , y are reals and a, b 1 c, d integers and the mod 1 operation means taking 
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the fractional part in (0,1) of the real number. In this paper we restrict to the 
study of periodic orbits, whose initial point, (x(0),y(0)) = ( q/p,q'/p ), is ratio¬ 
nal, with common denominator p. These orbits are then equivalently generated 
by the transformation 


( x(t + l)\ 


fa b\ f x{t)\ 
\c d) \y(t) ) 


mod p , 


(3) 


where x, y are integers belonging to [0,p— 1] (the field Z p ), thanks to the mod p 
operation. 

Moreover, we restrict to prime values of p, which allow for the existence 
of orbits with maximal period p — 1. Let A = (TrA) 2 — 4, if an integer q = 
0,..., p— 1 exists such that q 2 = A mod p, then matrix A admits two eigenvalues 


A± = 2- 1 (Tr A±q) , 


(4) 


where 2 1 is the multiplicative inverse of 2 in the field Z p . From these properties 
we can directly derive the equations for linear modular eigenspaces (LME) 


y = n±xmod p , 


(5) 


where 

n± = —b~ 1 (a — A±) , (6) 

and (x,y) £ Z p x Z p . 

A ’primitive root’ of p, is a positive integer g for which the smallest value of 
s such that g s = 1 mod p is s = p— 1. It has been observed □ that ordered pairs 
( g n ,g n+l ) in Z p x Z p with varying n = 1 ,,p— 1 and the fixed parameter l 
chosen in the set {1,... ,p— 1} distribute “orderly” or “chaotically” on the lattice 
of points Z p x Z p depending on the chosen value of l. This result is relevant 
for the understanding of correlations in the generation of random numbers by 
modular transformations of integers. 

These ordered pairs can also be seen as successive points over an orbit of 
map (||) with initial value ( g k , g k+l ) {k £ {0,... ,p — 1}) belonging to the LME 
(j5j). Therefore, choosing the value of l and the primitive root g is equivalent to 
fixing the elements of matrix A in (j3|) and the initial condition, in such a way 
that the following consistency relation is obeyed 

bg l - g + a = 0 . (7) 

In this paper we propose to study the properties of the pairs, or equivalently 
of the orbits of (3|), as a function of the parameter n± which fixes the LME, 
and we establish a relation of ’similarity’ among these orbits after introducing 
discrete rotations of the point lattice (Section 3). This allows us to overcome 
a difficulty encountered in Ref. □, where to the slightest variation of parameter 
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l corresponded a drastic change in the behavior of the orbit, from “order” to 
“chaos” and back, with no apparent regularity in the change. Using our approach 
it is possible to classify the different behaviors of the orbits according to the 
equation of the LME once p and g are given. 


3. Linear modular eigenspaces 

The study of LME equations (see formula (jj|)) as the parameter n± is varied 
is useful to determine the graph of linear periodic orbits of a specific dynam¬ 
ical system fixed by the values of p and g. The value of n± varies in the set 
{0,1,. -. ,p — 1}. Eq. (|) can then be rewritten as 

y = n±x - hp , (8) 

with h = 0,1,..., n± — 1 and [hp/n±\ + 1 < x < [(h + 1 )p/n±\. This is the 
equation of an ensemble of n± parallel segments, whose slope is n±. If n± 
divides p — 1, then the number of points of the LME over each segment is the 
same. Periodic orbits over a LME contain all points of the LME, but visit them 
irregularly, jumping among the segments of the broken line (|j). 

Let us give an example. Consider the matrix 

(' !) ■ < 9 > 

whose eigenvalues are A± = 2 ± 2 ~ l q, with q 2 = 12. Let us consider the prime 
p = 13 and the primitive root g = 6. There are two values of q whose square is 
12 in the field Z 13 , i.e. q = 8 and q = 5; we choose here q = 8 and then n+ = 3, 
n- = 8, A+ = 6 and A_ = 11, the latter being also primitive roots. Taking the 
initial point on the LME, e.g. (a;(0), y(0)) = (1,3) one generates pairs which 
remain on LME, but jump among the 3 segments in which the LME breaks up 
(see Fig. 1); the orbit has period 12. With an orbit of the same period one can 
obtain a more homogeneous filling of the point lattice. 

The best one can obtain for these values of p and g is for n+ = 5. A possible 
choice of matrix A is obtained by the following procedure. First choose a value 
of q which respects the condition q 2 = (a + d ) 2 — 4; then, using equations 

2 (a + n±b) = a + d±q (10) 

ad — be = 1 , (11) 


determine two of the elements of A in terms of the others (e.g. the diagonal 
terms a and d). In our case one obtains the matrix 


3 

11 



( 12 ) 


whose eigenvalue A+ = 11 (A_ = 6) with q = 5 (the other choice of q with 
respect to the previous example). In this case, since the prime p = 13 = 3 2 + 2 2 
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Figure 1: LME (dotted line) and a periodic orbit on it for p = 13, g = A+ = 6, 
n+ = 3. 



is the sum of squares, the square S p = [0,p — 1] x [0,p — 1] is tiled by p squares 
of area p (see Fig. 2), which realizes the most uniform filling. 

This recipe can be generalized to larger values of p which are sums of squares; 
the square of the diagonal of the tiling square is 2 p and can be written in terms 

of n .|_, 

2p = n\ + 1 , (13) 

which gives n+ = \J2.p — 1. This is the optimal choice for uniform filling, but it 
is valid only for specific values of p. In the following, we will find a more general 
rule valid for any p (see formula (|l8|)). 

Let us now introduce a ’similarity’ relation among LMEs. We say that two 
LMEs are similar if one can obtain a LME from the other by rotating it around 
the point (p/2, p/2) (intended as a real number) on the square S p . Let us 
consider the LMEs 

y = rriiX mod p (14) 

The similar LMEs are at most four, with (mo,mi,m 2 ,m 3 ) such that momi = 
lmod p, m 2 = P — mo and m 3 = p — m\. The LME with slope to,; is obtained 
rotating by ?'7 t/ 2 the LME with slope too- The proof of this property is trivial. 

4. Uniformity probes:entropy and geometric index 

In order to characterize the distribution of the points of the periodic orbit on 
the square S p , and its variation with n±, we define the ’entropy’ of the periodic 
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Figure 2: Same as Fig. 1 but n+ = 5, which realizes a more uniform filling of 
S 13 = [0,12] x [0,12]. 

orbit as follows. Let us divide the square S p into square cells of side 6 . Let Nk 
be the number of points of the periodic orbit which fall inside the k— th cell and 
T p = p — 1 the period. The entropy is 

M 

Tj = - Wk Inwk , (15) 

k =1 

where M = {p/5) 2 is the total number of cells and Wk = Nk/T p is the fraction 
of points of the orbit inside the k— th cell. To avoid the presence of empty cells 
(cells without points of the orbit) in the homogeneous case, 5 must be chosen 
such that M = {p/5) 2 = T p = p — 1, then 5 = p/yjp — 1. As a consequence 

0 < r? < ln(p — 1) . (16) 

Another simple probe of homogeneity is the ratio of the sides, L max and L m i n , 
of the elementary cell of the lattice of points formed by the orbit in S p , 

r _ Lmaa ' ( 17) 

Lmin 

Then r > 1 and the r = 1 case corresponds to the most uniform square elemen¬ 
tary cell (which is reached only for the primes which are sums of squares). The 
supremum value of r is r sup = p/\J 2, which is attained when the points of the 
orbit are all on a single segment (it is computed by considering the nearest S p 
square). 
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In Fig. 3 we show the dependence of rj on n, the slope of the LME. The chosen 
prime is p = 1999 and the primitive root is g = 3. The minimal entropy value, 
which is reached when the points of the orbit are placed on a single straight 
segment, is easily computed to be rjmin = 1/2 ln(p — 1) (equal to 3.799... in 
this case) and is obtained when n = 1,1998 (n = 1 ,p — 1 in general). The 
maximal entropy value rj ma x = 7.599... is reached for n = 45,844,1155,1954; 
these values are related through the similarity relation introduced in Section 3. 

It is remarkable that the homogeneity ratio r (see eq. [Tt]) gives a perfectly 
equivalent information, as shown in Fig. 4. We have verified that the agreement 
between the two homogeneity probes gets better as the value of the prime p 
increases. 
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Figure 3: Entropy of periodic orbits with p = 1999 as a function of n, also shown 
with dotted lines are the maximal, 7.599 ..., and minimal 3.799 ... entropies. 



Which is the optimal n opt value corresponding to the highest homogeneity? 
A simple heuristic argument suggests that it must be of the order of the side of 
the elementary square 6 = p/y/p — 1. In fact, numerical experiments suggest, 
with high degree of confidence, the following conjecture 


Tiopt — maxint( 


Vp - 1 


)■ 


(18) 


Fig. 5 shows the agreement of this conjectured value with the one determined 
numerically. Of course, the similarity conditions relate this value with three 
others. 


5. Random number generation 

One can use the transformation defined in (|^) by the LME to generate recur¬ 
sively points at random with coordinates ( X , Y) in the square S = (0,1] x (0,1] 
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Figure 4: Geometric uniformity factor r for p = 1999 as a function of n. The 
supremum r sup = 1413.506... is reached for n = 1 and for n = 1998 when the 
orbit collapses onto a segment. The supremum and the minimal r = 1 values 
are shown by the dotted lines. 


as follows 


X{t) = 

X{0)S ‘ mod 1 


P 


Y(t ) = 

S ' XW mod 1 . 

(19) 


once one has chosen the primitive root g and the value l such that g l = n±. 
X(0) £ Z p is the seed of the generator and can always be written in terms of 
the primitive root X(0) = g to with to £ Z p . 

In the previous section we have found the value of n± ([jj]), and consequently 
of l, which, for fixed g, gives the most uniform distribution of points on the 
square S if the full orbit with initial value (X(0) = g to , Y(0)) is generated. We 
want now to test whether the recursion © can also be used as an efficient 
generator of random pairs (X(t),Y(t)) with the appropriate choices of g and 
l. To test the efficiency we have used a method recently developed in Ref. 0 to 
characterize ergodicity of tranformations of the plane into itself. These authors 
introduce a random model where the square S is divided into N cells and then 
points are thrown into the cells each at a time with no correlations. At each 
step t of the process one has equal probability it = IV -1 , of choosing any cell. 
One can then define the probability distribution Pt(k) that k cells are occupied 
at step t. The first two moments of Pt{k) are computed in 0 


N 

p{t) = ^(/C7r)P t (fc) = 1- (1-71-)* 

k =1 


( 20 ) 
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Figure 5: Conjectured (line) and numerically determined (points) n op t as a 
function of the prime p. 



N 


S(t) = ^(fc7r) 2 P t (fc) = 1 — (2 — 7r)(l — 7r) 4 + (1 — 7r)(l — 27 t)* . (21) 


fc=l 


The first moment p[t) is simply the density of occupied cells at step n and can 
be approximated for large N as 


p{t) « 1 - exp - 


t 

N 


( 22 ) 


The second moment gives the standard deviation a(t) = y/S(t) — p 2 {t) at step 

t. 

In Fig. 6 we plot the numerically computed deviation with respect to the the¬ 
oretical result of the random model (|2^) for the random number generator ( |l9| ) 
with p = 2 31 — 1, g = 7 5 and the optimal value n op t (18) for the slope, together 
with the deviation for a commonly used random number generator, ranO □. 
The ±<r(t) values are also reported and show that generator (jid|) systematically 
overestimates the theoretical curve, while ranO fluctuates statistically around 
the theoretical value remaining within the ±<x(t) region. 

Our generator dH) has therefore the tendency to produce points on the 
square S which overfill different cells with respect to the random filling. This 
must be an effect of the temporal correlations which are present in ^) and 
tend to distribute points on parallel straight lines. Therefore, we have modified 
( |l9| ) imposing a variation at each step of the slope of these lines. Hence, the 
generator is changed as follows 


X(t) = 


*(o )g t 


mod 1 
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Figure 6: Deviations from the theoretical curve of the density of occupied 
cells in the random model for random series generated with map (jig) (full line) 
and ranO (dotted line). The ±cr(f) • 10 3 lines are also reported (broken lines). 


N = 10 6 . 


r(f, = d^mocll (23) 

= r!ifflo modl , 

P 

with the seed (X(0),r(0)) £ Z p x Z p . As previously, the seed can always be 
written in terms of the primitive root (X(0) = g to ,r( 0) = g s °) with to and so 
belonging to Z p . 

The randomization of the slope r reduces the correlations and now the gen¬ 
erated random series falls within the error curves ±cr(t) (see Fig. 7), showing 
the same quality of ranO. 

Although we do not claim here that generator (^3[) might be competitive 
with others, more widely used, random number generators, we are convinced to 
have revealed an important feature of modular generators, i.e. cell overfilling 
tendency. The recipe we have proposed to balance this effect is not be unique 
and other solutions are certainly possible. 

6. Conclusions 

We have studied the geometrical properties of the orbits of modular trans¬ 
formations on integers. We have solved a puzzle proposed in Ref. □ related to an 
“order” to “chaos” transition as a parameter of the transformation is varied by 
showing that indeed the true control parameter is the slope of linear modular 
eigenspaces; we have then determined which of the values of this latter control 
parameter give an “ordered” (i.e. distributed over a few segments) or “chaotic” 
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Figure 7: Same as Fig. 6 where (Jh|) has been substituted by (23). 


(i.e. spread over the lattice) orbit by looking at uniformity probes. 

Particular attention has been devoted to uniform (“chaotic”) orbits, which 
we have characterized using entropic and geometric parameters. Moreover, we 
have verified that, during the time evolution, the most “geometrically” uniform 
orbits do not correspond in general to the most “statistically” uniform, as defined 
through a specific model introduced in Ref. 0 . 

With the aim of generating statistically uniform orbits we have introduced 
a new generator of pseudo-random pairs, where the slope of the linear modular 
eigenspace is varied during time evolution. The recipe we propose to generate 
pairs ( X , Y) of pseudo-random real (32-bit floating point) numbers with uniform 
distribution in the square S = (0,1] x (0,1] is by following recursion (t = 0,1,...) 


m = 

X(0 )g t 

--— mod 1 

P 

r{t) = 

^ mod 1 


P 

Y(t) = 

r(t)X( ‘> mod 1 


(24) 


where p is a prime (e.g. the typical value p = 2 31 — 1), g a primitive root 
(e.g. g = 7 5 = 16807) and the seeds X(0),r(0) are integers from 1 to p — 1. 
The comparison of this generator with others commonly used is encouraging, as 
far as the uniformity properties are concerned. The method can also be easily 
extended to generate fc-tuples of random numbers in fc-dimensional spaces. 
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